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In the face of mounting numerical evidence, Metlitski and Grover arXiv: 1112.5166 have given 
compelling analytical arguments that systems with spontaneous broken continuous symmetry con¬ 
tain a sub-leading contribution to the entanglement entropy that diverges logarithmically with 
system size. They predict that the coefficient of this log is a universal quantity that depends on 
the number of Goldstone modes. In this paper, we confirm the presence of this log term through 
quantum Monte Garlo calculations of the second Renyi entropy on the spin 1/2 XY model. Devising 
an algorithm to facilitate convergence of entropy data at extremely low temperatures, we demon¬ 
strate that the single Goldstone mode in the ground state can be identified through the coefficient 
of the log term. Furthermore, our simulation accuracy allows us to obtain an additional geometric 
constant additive to the Renyi entropy, that matches a predicted fully-universal form obtained from 
a free bosonic field theory with no adjustable parameters. 


Introduction - In condensed matter, the entanglement 
entropy of a bipartition contains an incredible amount of 
information about the correlations in a system. In spa¬ 
tial dimensions d >2, quantum spins or bosons display 
an entanglement entropy that, to leading order, scales 
as the boundary of the bipartition m- Subleading to 
this “area-law” are various constants and - particularly 
in gapless phases - functions that depend non-trivially on 
length and energy scales. Some of these subleading terms 
are known to act as informatic “order parameters” which 
can detect non-trivial correlations, such as the topologi¬ 
cal entanglement entropy in a gapped spin liquid phase 
mz]. At a quantum critical point, subleading terms con¬ 
tain novel quantities that identify the universality class, 
and potentially can provide constraints on renormaliza¬ 
tion group flows to other nearby fixed points [HHIl]- 
In systems with a continuous broken symmetry, ev¬ 
idence is mounting that the entanglement entropy be¬ 
tween two subsystems with a smooth spatial bipartition 
contains a term, subleading to the area law, that diverges 
logarithmically with the subsystem size. First observed 
in spin wave |15j and finite-size lattice numerics |16] , the 
apparently anomalous logarithm had no rigorous expla¬ 
nation until 2011, when Metlitski and Grover developed 
a comprehensive theory m- They argued that, for a 
finite-size subsystem with length scale L, the term is a 
manifestation of the two long-wavelength energy scales 
corresponding to the spin wave gap, and the “tower of 
states” arising from the restoration of symmetry in a fi¬ 
nite volume [HHH]. Remarkably, their theory not only 
explains the subleading logarithm, but predicts that the 
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FIG. 1. Schematic energy level structure of the low energy 
tower of states for finite-size systems with spontaneous break¬ 
ing of a continuous symmetry. The correction to the entangle¬ 
ment entropy may be approximated by the log of the number 
of quantum rotor states below the Goldstone gap, Qa, which 
is represented by the states within the dotted box. 



coefficient is directly proportional to the number of Gold¬ 
stone modes in the groundstate. Furthermore, describ¬ 
ing a Goldstone mode with a free scalar field theory al¬ 
lows them to predict the value for an additional additive 
geometric constant, which is fully universal and should 
therefore be the same across a wide range of continuum 
theories and lattice models. 

In this paper, we confirm these predictions in a striking 
way, through large-scale quantum Monte Garlo (QMG) 
simulations on a spin-1 /2 XY model on the square lattice. 
By employing an extended-ensemble generalization of a 
ratio method [22H24] . we are able to carefully converge 
the mutual information of the second Renyi entropy to 
its low-temperature value. There, finite-size scaling re¬ 
veals the coefficient of the subleading logarithmic term to 
precisely match the prediction of Metlitski and Grover, 
identifying the lone Goldstone mode in the theory. Re¬ 
markably, our simulations are accurate enough to also 
measure the universal additive geometric constant, which 













2 


is consistent with the predicted relationship to the con¬ 
stant which appears in a lattice-regularized free scalar 
field theory m 

Entanglement entropy in the tower of states - To ob¬ 
tain a qualitative understanding of the origin of the log¬ 
arithmic correction in Metlitski and Grover’s theory, it 
is simplest to first envision decoupling the two spatial 
subsystems, A and B , that define the entangled biparti¬ 
tion. The low-energy degrees of freedom in each subsys¬ 
tem can be described by an 0{N) rotor (TV = 2 for the 
XY model), representing the direction of the order pa¬ 
rameter. Here we are only allowing global fluctuations of 
the order parameter within each subsystem, such that we 
may approximate the state of A and B each as a single in¬ 
dependent quantum rotor. The effective Hamiltonian of 
each subsystem is H = f 21, where lA is the total an¬ 

gular momentum operator with eigenvalues £{£+!) and 
/ is the effective moment of inertia which is extensive, 
proportional to the magnetic susceptibility y: / ~ yL'’* 
in d spatial dimensions [5^ . Thus the energy scale of the 
tower Atow = 1/x^^ vanishes with the system volume, 
faster than any other energy scale. The eigenstates of 
result in the famous “tower of states” observed routinely 
in computational studies of systems with continuous sym¬ 
metry breaking in a finite volume [IMH- 

The interaction between A and B which aligns the sub¬ 
system order parameters may be introduced via a Gold- 
stone mode Hamiltonian Hq which couples the two ro¬ 
tors. The energy scale of the Hq is the Goldstone mode 
gap Aq which is the scale of the lowest energy spin 
waves. Since Aq ^ cjL where c is the spin-wave velocity, 
Aq ^ Atow in the thermodynamic limit for d > 1 . In 
the limit Aq —>■ oo, there are no relative fluctuations in 
the order parameter between subsystems, and A and B 
act as a single rigid rotor. For finite Aq, there will be 
relative fluctuations between the subsystems order pa¬ 
rameters due to the zero point fluctuations of Hq. 

To estimate the entanglement entropy contribution 
from the tower of states, we can count the number of 
“accessible” states of subsystem A, when the total 
system is in the ground state, and use S'tow logH^. 
In the limit Aq —>■ cx) and the rotors are rigidly cou¬ 
pled, the ground state is the ground state of the total 
system tower of states with zero total angular momen¬ 
tum: £ab = 0. In this case all states in the A subsystem 
tower are accessible to the ground state, as each state in 
A can be paired with an appropriate state in B to form 
a state with nonzero overlap with the £ab = 0 state. 
However, as discussed above, by including a finite Aq 
and thus allowing relative fluctuations of the subsystem 
order parameter between A and B, the fluctuations in 
the subsystem angular momentum are finite and deter¬ 
mined by the ratio of the energy scales: (L^) ^ Aq/Atow 
m In fact, the reduced matrix of the subsystem takes 
the form of a thermal density matrix with an effective 
“entanglement Hamiltonian” given by Htow and the “en¬ 


tanglement temperature” given by Ag m-: the resulting 
tower of states structure in the entanglement spectrum 
has been seen in numerics |26H29) . Thus the inclusion 
of Goldstone modes cuts off the accessible states of the 
subsystem to those with an energy below the spin wave 
gap, as illustrated in Fig. 

As an example of this mechanism, consider the case 
of TV = 2 (valid for our XY model simulations below). 
Here the rotors have a single component of angular mo¬ 
mentum £^ and the orientation of the rotors is described 
by a single angle 9. For Aq —)■ oo the ground state has 
^AB = Oi which has nonzero overlap with states of equal 
and opposite £^ in each subsystem: \£\ = £,£g = —£); 
consequently all |£a) state are accessible in this limit. We 
may include the effect of the lowest Goldstone mode by 
treating the dynamics of the relative angle between sub¬ 
systems 9s as a single harmonic oscillator with frequency 
Ag and moment inertia Is with an effective 

Hamiltonian 

Hg = + \l5Al9l ( 1 ) 

Here, the fluctuations in the relative angular momen¬ 
tum Ls are given by the ground state fluctuations of a 
harmonic oscillator: {Lf) ~ Is Aq/2 ^ Ac/Atow The 
key point here is that because the order parameter is 
canonically conjugate to the rotor angular momentum, 
increasing the relative fluctuations in the order parame¬ 
ter reduces the fluctuations in L^. Thus, allowing relative 
fluctuations of the order parameter between subsystems 
effectively cuts off subsystem rotor states that are ac¬ 
cessed in the ground state at order £ ~ {Aq/A towY^"^ - 
a relationship that holds for all TV m- 

We may therefore estimate VIa by counting the num¬ 
ber of states (in A’s tower of states) that lie below Ag. 
For systems with 0{N) symmetry, the tower of states 
is described by a rotor living on an Nq = TV — 1 di¬ 
mensional sphere, where Nq is the number of Goldstone 
modes. The dengeneracy of each energy level is of order 
£Ng-i^ We then may estimate the total number of states 
below Aq by integrating the degeneracy up to the cutoff 
4o = (AG/Atow)'/": 

rEo /A \ Ng/2 

Ha^ d££^--^^(^] . (2) 

Jo V ^tow / 

Using the relation y = Ps/c^ from hydrodynamic spin- 
wave theory where ps the stiffness m, the entanglement 
entropy correction due to the tower of states becomes 

5tow~^log(^T"-'). (3) 

We see that the logarithmic correction to the area law 
arises due to the quasi-degeneracy of accessible bulk sub¬ 
system states, that scales as a power law in L for systems 
with spontaneously broken continuous symmetries. This 
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contrasts with the leading-order area law, arising from 
the exponential scaling of the number of local boundary 
states with the boundary area. Clearly, the prefactor 
of the logarithmic correction is a universal number that 
simply counts the number of Goldstone modes. 

Quantum Monte Carlo Procedure - In order to ex¬ 
amine the effects of broken continuous symmetry on en¬ 
tanglement, we implement a highly-efhcient Stochastic 
Series Expansion (SSE) QMC algorithm [3TVI33] for the 
2D spin-1/2 XY model, H = + 5f5f). 

This model is known to realize a ground state where 
the U(l) symmetry is spontaneously broken, resulting 
in one Goldstone mode. The finite-temperature SSE al¬ 
gorithm uses a version of the directed-loop updates spe¬ 
cialized to the XY model [34]. To measure the entan¬ 
glement entropy, we employ a replicated simulation cell 
Z[A,2,T], which gives access to second R&yi entropy 
Sn = log [Tr/9^]/(l — n), with n = 2. This is done 
through the ratio Trp^ = Z[A, 2, T]/Z[A=0, 2, T], where 
Z[A=0,2,T] = Z[T]'^] the square of the unmodified par¬ 
tition function [35]. One therefore needs to evaluate a 
difference in free energies or, alternatively, a ratio of par¬ 
tition functions, between the replicated and unmodified 
simulation cells with QMC. Several options are avail¬ 
able, such as thermodynamic integration [^ and Wang- 
Landau sampling [36]. Alternatively, methods such as 
the ratio method allow one to compute the par¬ 

tition function ratio directly. As described in the next 
section, for this model, careful convergence to low tem¬ 
perature is required. Hence, in the Supplemental Ma¬ 
terial, we develop a highly-efhcient variant of the ratio 
method, dubbed the extended ensemble ratio method, for 
the XY model. 

Results - We perform large-scale calculations to ob¬ 
tain the second Renyi entropy of the XY model on Lx L 
square lattices with periodic boundary conditions. The 
resulting torus is partitioned into two cylindrical regions 
of linear dimension L x and L x [L — L^). For sys¬ 
tem sizes L = 8, 12,16, separate Hnite-temperature tests 
are performed to explore the convergence of S 2 - In agree¬ 
ment with the expected scaling of the tower-of-states gap, 
those results conhrm that the convergence temperature 
scales approximately as 1/L^. We use this to estimate 
the convergence temperatures for L = 20, 24, 28, 32. 

Since SSE QMG simulations necessarily run at hnite 
temperature, very small thermal contributions to S 2 are 
expected, which we observe to signihcantly affect our 
hnite-size scaling analysis below. Fortunately, this ther¬ 
mal contribution can be essentially eliminated by em¬ 
ploying the mutual information, l 2 {A : B) = S' 2 (A) -f 
82 ( 3 ) — S 2 {AB) [35]. Assuming that bulk (volume-law) 
contributions from the two subsystems A and B approx¬ 
imately cancel the bulk term S 2 {AB) in I 2 , the scaling 
form described by Metlitski and Grover becomes: 

(4) 



FIG. 2. (color online). The mutual information of the S' = 1/2 
XY model as a function of /3 = J/T. Solid lines are obtained 
through thermodynamic integration from /3 = 0, with statisti¬ 
cal errors estimated by the shading. Square points with error 
bars are data obtained at a fixed /? using the extended ensem¬ 
ble ratio method, described in the Supplemental Material. 


Here, a is a non-universal constant, ps and c are the 
spin-wave stiffness and velocity, and 7ord is the geometric 
constant that depends on the aspect ratios of the cylin¬ 
ders A and B m- Note, since non-universal (cutoff) 
dependences are all contained within ps and c, this geo¬ 
metric constant remains fully universal. For the spin-1/2 
XY model on the square lattice, ps = 0.26974(5) J and 
c = 1.1347(2) J were obtained from Ref. [37]. 

Figure [^ illustrates a representative convergence test 
for different system sizes. The mutual information peaks 
at temperatures above the Kosterlitz-Thouless transi¬ 
tion of {T/J)kt = 0.343 (which can be detected by 
the crossing of the finite-size curves; see Ref. [35] )• For 
T/J< (T/J)kt, the mutual information reaches a min¬ 
imum (at J/T = /3 « 4 in Fig. before undergoing 
a slow rise. This rise continues until the approximate 
ground state is reached, for temperature below the finite- 
size scaling gap, which for system sizes larger than L = 8 
occurs for /3 > 100. Thus, although the method of ther¬ 
modynamic integration is useful to produce the general 
shape of the I 2 curve for a wide range of temperatures, 
it is difficult to control the systematic error introduced 
by numerical integration at low temperatures for L > 12. 
Therefore, data used in the below fits was converged at 
very low temperatures using the extended ensemble ratio 
method, described in the Supplemental Material. 

Figure [^ illustrates the resulting temperature- 
converged mutual information for a variety of system 
sizes, as a function of the “width” of the cylindrical re¬ 
gion, Since, for a subsystem A and its complement 
B, Sa = Sb only at T = 0, the symmetry of the en¬ 
tanglement entropy about L^/L =1/2 provides a sensi¬ 
tive test of temperature convergence. Employment of the 
“bare” Renyi entropy results in a very slight asymmetry 


l 2 = aL + No log{Lps/c) + 27ord- 
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FIG. 3. (color online). The mutual information as a function 
of torus aspect ratio, for the lowest temperatures examined 
for each system size. The corresponding /3 are 184, 368, 736, 
1150, 1650, 2300, 3200 ordered from the smallest to the largest 
system size. Vertical dashed lines are the aspect ratio values 
employed in the fitting in Fig. 


in the curve; use of I 2 restores this symmetry producing 
high-quality data that can be fitted using Eq. (|^. 

The results of this analysis are illustrated in Fig 
Here, I 2 is calculated at various aspect ratios (the verti¬ 
cal cuts in Fig. and fit to the functional form Eq. 
Specihcally, to extract the coefficient of the sublead¬ 
ing logarithm, the mutual information was fit to I 2 = 
aL + b\og{Lps/c) -\- d, where a, b and d are adjustable fit 
parameters. As illustrated in Fig. (a), there is defina- 
tive evidence for the existence of a logarithm; further¬ 
more, independent fits for the four aspect ratios studied 
each give Nq = 1 to within error bars. 

Even more striking, we are able to extract the uni¬ 
versal shape-dependence of the geometric constant 7ord- 
To do so, fits were performed to the functional form 
I 2 = aL + \og{Lps/c) + 27ord + d/L, where Nq is fixed at 
unity in order to remove one parameter from the analy¬ 
sis. Thus calculated, 7ord for TV = 2 in two dimensions 
can be compared via a zero-parameter fit to the sublead¬ 
ing constant term 7free calculated in a free scalar field 
theory m through the relation 7ord = 7free + | log(27r), 
valid for the second Renyi entropy. The free field result 
yfree, which depends on the aspect ratio can be 

calculated numerically for free bosons on the lattice us¬ 
ing the correlation matrix technique (as in Ref. [39]). As 
illustrated in Fig. |^(b), the resulting theoretical curve is 
in excellent agreement with our QMC results for 7ord- 

Discussion - In this paper, we have employed large- 
scale quantum Monte Carlo (QMC) simulations on the 
spin-1/2 XY model to demonstrate the presence of a log¬ 
arithmic correction to the Renyi entropy due to sponta¬ 
neous breaking of continuous symmetry. This term arises 
from the presence of two infra-red energy scales in the 
problem; the spin-wave gap, and the “tower of states”. 




FIG. 4. (color online). The top figure displays a three pa¬ 
rameters fit to the functional form aL -\- b\og{Lps/c) + d for 
different torus aspect ratios, where b gives a value for Nq. An¬ 
other three parameters fit to aL + \og{Lps/c) -I- 27ord -f d/L 
(not shown) is employed to extract the geometrical constant 
7ord. The resulting Nq and 7ord are shown alongside the 
theoretically-predicted values at bottom. 


The coefficient of this logarithm is predicted in Ref. m 
to be NQ{d — l)/2, where Nq is the number of Gold- 
stone modes and d the spatial dimension. We confirm 
this prediction in a striking manner, through finite-size 
scaling studies on the square lattice XY model, recover¬ 
ing to high precision the expected Nq = 1. In order to 
do so, a QMC algorithm was developed to sample the 
second Renyi entropy at a fixed temperature with high 
efficiency (described in the Supplemental Material). 

Remarkably, in addition to confirming Nq = 1, our 
simulation technique is able to converge the value of an 
additional additive geometric constant 7ord, which is fully 
universal since all short-distance physics is confined to 
the (known) spin wave stiffness and velocity, contained 
within the argument of the logarithm. The resulting 7ord 
has a functional dependence on the geometric aspect ra¬ 
tio of the entangled bipartition. This function matches, 
to within error bars, that calculated using a free scalar 
field theory regularized on a toroidal square lattice, with 
no adjustable parameters. 

This is a rare and striking example of complete quanti¬ 
tative agreement between a universal quantity calculated 
in continuum field theory and finite-size lattice simula¬ 
tions. The ability of the Renyi entanglement entropy to 
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mediate between continuum theory and lattice simulation 
illustrates the utility of such geometric quantities as an 
alternative paradigm to study correlations in condensed- 
matter systems. In this case, with the full understanding 
of the universal structure of the entanglement entropy in 
the presence of a spontaneously broken continuous sym¬ 
metry, the door is now open to the direct search and de¬ 
tection of Goldstone modes in a large variety of systems, 
through the widely accessible Renyi entropies. 
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Supplementary material for “Detecting 
Goldstone Modes with Entanglement 
Entropy” 


defined with two arbitrary regions A and A'; namely 
Z[A',2,T]/Z[A,2,T]. As with all ratio methods, the 
QMC estimator in this case is based on a simple iden¬ 
tity [S4] : 


PRELIMINARIES 


This Supplemental Material describes a quantum 
Monte Carlo (QMC) algorithm for estimating the Renyi 
entropy, 

^n = 7^1og[Trpl], (SI) 

1 — n ■' 


using an extended-ensemble version of the so-called “ra¬ 
tio method” first applied in this context by Humeniuk 
and Roscilde m- In the present work, the algorithm will 
be specialized to the case of the second Renyi entropy 
(n = 2), although extensions to other n are relatively 
straightforward. Similarly, although we concentrate on 
applications to the spin-1/2 XY model, generalizations 
to Hamiltonians with other symmetries are possible. 

QMC estimators for the second Renyi entropy S 2 rely 
on the fact that the trace over powers of the reduced den¬ 
sity matrix can be related to a ratio of partition functions: 


Trpi 


Z[A,2,T] 

Z[A=0,2,Ty 


(S2) 


where the numerator is the partition function of a multi- 
sheeted Riemann surface [S2] , (a “replicated” QMC sim¬ 
ulation cell) and the denominator is the square of the 
regular partition function, Z[A=0,2,T] = Z[Ty. In the 
replicated case, the geometry of the entangled region A 
dictates the “boundaries” of the d -I-1 dimensional QMC 
simulation cell; world-lines in region A are periodic in 
imaginary time with period 2/3 (or n/3), while in region 
B, 2 (or n) independent replicas exist with periodicity /3. 

Since the logarithm in Eq. (SI) reduces the calculation 
of S 2 to the difference in free energies between systems 
described by the two partition functions, thermodynamic 
integration or Wang-Landau techniques can be used to 
devise QMC estimators. In contrast, ratio methods (not 
to be confused with the related “ratio trick” coined in 
Ref. [S3]) give an estimator for the ratio of partition func¬ 
tions, and act as valuable alternatives to explicit calcu¬ 
lation of free energies. The most distinct advantage of 
ratio methods is their ability to calculate 5'„ directly at 
a given fixed temperature. They can however be ineffi¬ 
cient for large entropies (large subregions A), and may 
require several separate simulations of different subregion 
geometries (that can be combined with the ratio trick) 
to combat this. 


EXTENDED-ENSEMBLE RATIO METHOD 


We begin by generalizing Eq. (S2) to the problem 


of calculation the ratio of replicated partition functions 


Za \ Wa(c) ]^ ’ 


(S3) 


where we have simplified our previous notation Za = 
Z\A,2,T\. The expectation value is taken in the Za en¬ 
semble and Wa'{c) and Wa{c) are the weights of a con¬ 
figuration c in the Za' and Za ensembles respectfully. 

The difficulty of applying this identity in a straightfor¬ 
ward manner is that it is valid only when the configura¬ 
tion space of Za', ^A', is contained within the configu¬ 
ration space of Za^ a ~ that is W{c) ^ 0 for any c. Un¬ 
fortunately, this condition is not automatically satisfied 
when the standard labelling variables are used in conven¬ 
tional QMC techniques. However, for the XY model, by 
introducing a new label, it is possible to recast the config¬ 
uration spaces of both ensembles such that Qa = ^a' ■ In 
other words, a common label c can be found that enumer¬ 
ates all possible configurations within both ensembles. 

Here we will consider a stochastic series expansion 
(SSE) representation of the partition function [SStiSTj . 
A term in the SSE expansion of Z can be represented 
by a spin state and an operator list, where the opera¬ 
tors are terms in the Hamiltonian. Alternatively, we can 
represent the same configuration by a linked list of ver¬ 
tices |S8]; an example of a set of vertices is displayed 
in Fig. |S1| We can therefore label the configuration by 
the list of links, I, and the set of vertex types v"^{l) for 
a given boundary condition defined by A as shown in 
Fig. |S2[ Here we have defined I to include all “internal” 
links between vertices and we link the “exterior” vertices 
to the boundaries of the operator string; we specifically 
do not close the links at the boundaries, as one would do 
in a regular d -|- 1 dimensional QMC simulation cell |S8] 

These two labels are sufficient to identify the configu¬ 
ration space of many models. Thus, a partition function 
can be expressed as the following double sum: 

^ W(u^(0), (S4) 

I V^{1) 

where W{v^{l)) is the weight of a configuration labelled 
by I and v^(l). 

This general expression can be simplified for the XY 
model. The model’s SSE vertices are displayed in Fig.|ST| 
With an appropriate choice of adjustable SSE constants, 
the weight of each vertex becomes equal to 1/2J [S8|. 
Hence, the weight of a vertex configuration does not 
depend on a particular combination of vertices in the 
list and, therefore, is completely defined by the length 
of the corresponding operator list alone. We get that 
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FIG. SI. (color online). The six types of vertices present in 
a SSE simulation of the XY model |S8| . The colored circles 
represent two possible states of a spin one-half. The horizontal 
solid bar depicts a two sites operator. The lines with arrow 
tips show all possible non-bounce moves for a given vertex 
when the entrance leg is the left-most bottom one. There are 
two of those moves for each vertex. 

Wxviv^il)) = Wxy{1) and Eq. ( [Sd] ) is simplified to: 

zr = (S5) 

I V^{1) 

The second sum counts the degeneracy of vertex con¬ 
figurations compatible with the boundary conditions be¬ 
tween replicas. Motivated by the search for a variable 
independent of those boundary conditions, we introduce 
a new label, s{v"^{l)), that enumerates all possible parti¬ 
tions of a vertex configuration labelled by (/, v"^{l)) into a 
set of non-overlapping segments. The following algorithm 
is used to construct a single instance of those segments: 

1. Pick an unmarked leg located on a boundary slice. 
Mark it as visited. 

2. By following the linked list, switch to a leg con¬ 
nected to it. 

3. The new leg belongs to a vertex. 

• If this vertex is unmarked, pick with an equal 
probability one of two possible non-bounce 
moves for this vertex and switch to the cor¬ 
responding leg. Mark this vertex as visited 
and store the move type. 

• If the vertex is marked, switch to the next leg 
by performing a move of the same type that 
was done before. 

4. Repeat steps 2-4 until a leg on a boundary slice is 
reached. 

By repeating this algorithm for all legs located on the 
boundary slices, all open segments are traced out. How¬ 
ever, it is possible that some of the legs located on the 
inner slices have been left unmarked after this procedure. 
In order to partition those remaining legs too, the closed 
segments (loops) need to be traced. This is achieved by 
adjusting two steps of the algorithm. Now in the first 
step, the choice of legs to be picked is extended to all 
interior legs. Once the initial leg is picked, the algorithm 
proceeds in the same way until it reaches the same leg 
again. Hence, the condition to terminate the execution 



FIG. S2. (color online). A three step conversion process of a 
vertex configuration v^{l) into (1) that preserves its seg¬ 
ment partition. The details of this process are illustrated in 
the text. At each step, a simulation cell composed of two 
replicas (top and bottom) is shown. The dashed vertical lines 
between adjacent legs that are located on different vertices 
form a linked list. The label I identifies different sets of those 
lines. The two types of small arrows placed next to the repli¬ 
cas’ boundary slices mark the boundary conditions along the 
dimension of the expansion: within a column, spins decorated 
with the same kind of arrows are connected. In this way, the 
first simulation cell’s region A is empty while the other two 
simulation cells’ region A! contains all three spins. Open col¬ 
ored solid lines trace out a segment partition of the first and 
the third simulation cells (7 segments total). Note that within 
the same cells, there is also a single closed segment, an inner 
loop, composed of four legs. In the second simulation cell, the 
open segments are merged by boundary connections to form 
a single cross-replica loop identified by the same color. Mis¬ 
matching boundary spins along this loop are flipped according 
to the algorithm presented in the text, resulting into a vertex 
configuration compatible with the new boundary conditions 
as displayed in the third simulation cell. 

of the forth step has to be modified appropriately. By 
construction, any two segments built in such a manner 
can never pass through the same leg and, therefore, are 
non-intersecting. 

The loops tracing continues until all legs are marked. 
By the end of this procedure every leg belongs to one sin¬ 
gle segment (closed or open). This constitutes a single 
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instance of the partitioning of a vertex configuration into 
a set of non-overlapping segments. An example of such 
partition is shown in the first cell of Fig. |S2[ Note that 
at each vertex, there are two choices how to proceed with 
the construction of a segment. Each of them leads to a 
different partition. [? ] Therefore, a simulation cell that 
contains Ny vertices can be partitioned in 2 ^'’ distinct 
ways. This is the range of the newly introduced label 
s. Since Ny just counts the number of vertices without 
discerning their types, the number of partitions for a par¬ 
ticular v"^{l) is determined by the I label only. This fact 


allows to rewrite Eq. (S51 in a new form: 




= E '■ (») 


v^{l) s{v-^{l)) 


This expression is obtained from Eq. (S51 by introducing 


a third sum via the substitution: 1 = 


where the new sum is performed over all different parti¬ 
tions of v^{l). 

It can be also shown that for any vertex configuration 
in A, v^{l), partitioned as s{v^{l)), there exists a vertex 
configuration in A', v'^ (1), with exactly the same par¬ 
titioning, that is s(u'^ (/)) = s{v"^{l)). The proof is by 


construction. If v^{l) and {1) were the same, the task 
is trivial. Otherwise, v"^{l) has to be modified in order 
to satisfy the boundary conditions of A!. An example of 


such process is displayed in Eig. Here, the first and 
thirds simulation cells represent v"^{l) and v'^ (1) corre¬ 
spondingly. The second cell depicts an intermediate step 
of the correctional procedure. Here, the open segments 
are connected into a loop along which the boundary spins 
mismatches are fixed one-by-one. Eurther details of the 
algorithm are given below. 

Proceeding column by column, consider each pair of 
boundary legs, (sj, s®) Ao be matched with respect to the 
new boundary conditions A'. If the legs align, 5° = s®, 
proceed to the next pair. Otherwise, randomly choose 
one of the two legs in the pair. Say it is s®. Since this 
leg is located on a boundary slice, it belongs to an open 
segment. Flip all legs belonging to this segment. Now, 
the original pair of legs is properly aligned, however there 
might be another mismatch at the other end of the seg¬ 
ment. Call the new pair (s}, S 2 ) where S 2 is the leg that 
has just been flipped as part of the open segment. By 
the same logic as before, sj must belong to an open seg¬ 
ment whose other end is identified as another boundary 
leg S 2 - if Si 9 ^ S 2 , flip this segment in order to align 
the (si,S 2 ) pair and move on to the next pair (s 2 ,sf). 
Otherwise, proceed to the same pair without flipping the 
segment. In this way, one-by-one pairs of boundary legs 
are aligned with respect to A boundary condition along a 
loop of open segments. An important subtlety occurs at 
the last step of this algorithm when the last pair (s^, S 2 ) 
is considered. Unlike previously, sj cannot be flipped if 


those legs do not align. An attempt to do so would en¬ 
tail another iteration of corrections with the same result, 
thus, initiating the algorithm in an infinite loop. 

However, this does not occur in the XY-model due 
to the special properties of its vertices. Notice that the 
only vertex move that connects two anti-aligned legs is 
the “switch-and-reverse” move (FiglSgi; this is the only 
move that reverses the vertical directionality of propaga¬ 
tion of the segment’s head. Consequently, once a segment 
tracing is initiated with the choice of a leg and its state, 
the spin state of the leg at the segment’s head is deter¬ 
mined by the vertical direction that the segment passes 
through the leg. On the last boundary connection, the 
direction of motion along the segment must be the same 
as the initial direction, and therefore the initial spin state 
at the head of the segment under construction is always 
the same its final state. We see then that for any seg¬ 
ment partition of v^{l), it is always possible to construct 
a v' (1) with the same segment partition. 

Now that we have shown the equivalence between any 
two configuration spaces constraint by boundary condi¬ 
tions A and A in terms of the segment partitions, we 
have achieved our initial goal to find a label s that can 
be used to apply the Eq. However, we still need to 
compute the weights of and (s). Upon a close 

inspection of the inner double sum in Eq. |S6[ we realize 
that the same segment partitions can be generated from 
different vertex configurations. This degeneracy can be 
exploited in order to replace the double sum with a single 
sum: 


(S7) 

I s{l) 

where the inner sum iterates over all unique segment par¬ 
titions for a given linked list and deg"^(s(^)) is the degen¬ 
eracy of the segment partition labelled as s(Z). As will 
be shown, the degeneracy depends on the boundary con¬ 
ditions between replicas and, thus, the superscript must 
be included. 

In order to calculate the degeneracy of a segment par¬ 
tition s(Z), connect open segments using boundary con¬ 
ditions A to form N^{s{l)) loops that cross boundary 
slices. In addition to those loops, there are also Ni{s{l)) 
inner loops (closed segments). Since all those loops do 
not intersect with each other by construction, the spins 
within them can be flipped independently. Each combi¬ 
nation of the loops’ flips leads to another vertex config¬ 
uration. Equivalently, all of those vertex configurations 
generate the same segment partition s(l)). There are in 
total N^{s{l))+Ni{s{l)) loops with each one being in one 
of two states: flipped or not-flipped. Each combination 
of those states corresponds to a different vertex configu¬ 
ration. In total, there are Qf such com¬ 

binations which constitutes the degeneracy deg"^(s(/)). 

With the last piece of the puzzle in our hands, we infer 
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the segment partition weight from Eq. (S7): 


Since the inner-loops are unaffected by the boundary 
conditions, upon the substitution of this weight in the 


Eq. (S3) their number drops out and an elegant expres¬ 


sion is obtained: 

Za 


= {2 


NC 


A- 


(S9) 


In practice, this estimator can be implemented in fewer 
steps that were required to prove its validity. It requires 
two routines. One routine traces out a random single 
segment partition s(l) for a vertex configuration v^(l) as 
was outlined before. In order to speed up the execution, 
it is not necessary to identify the closed segments. The 
end product of this routine is to associate the pairs of 
boundary spins that are connected via open segments. 
Once this step is done, the second routine takes the set 
of those pairs together with replicas’ boundary conditions 
as its inputs. Its task is to count This routine is 
executed for both A and A! with the same open segments. 
In the end, N^'{s{l)) and N^{s{l)) are known and the 
estimator can be evaluated according to Eq. (S9). 


BENCHMARKS 



To illustrate the efficiency of the new estimator, its raw 
measurements are compared to the original ratio method 
m in Fig. |S3[ on the 2D XY model of interest in the 
main text. Here, the deterioration of both estimators’ 
statistics is seen as the difference between the partition 
functions’ regions A, AA=A' — A, grows large. As the 
reference values, we employ the results obtained from the 
ratio trick |S3| , which constitute a compilation of the ex¬ 
tended ensemble (EE) ratio method results from many 
different Monte Carlo simulations, each executed with 
AA=1. Note that the original ratio method results are 
based on five times more Monte Carlo sweeps that were 
involved to produce the EE ratio method results. Even 
with such advantage, the ratio method statistics becomes 
increasingly poor towards AA=22. After this threshold, 
the estimator seems no longer capable to capture any 
meaningful statistics within the running time of its sim¬ 
ulations. The performance of the EE is strikingly better. 
Even when AA=64, the largest possible increment in the 
system, it produces an accurate result with a precision 
comparable to the ratio method precision at AA=10. 
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FIG. S3. The comparison of ratio methods efficiencies in a 
8x8 system with periodic boundary conditions at /3 = 8. The 
ratio of partition functions measurement is plotted against the 
size difference between their corresponding regions A, AA. 
The ratio method completely fails for AA > 22, and thus the 
data is not shown on the plot. The values obtained from the 
ratio trick serve as a reference. The statistical error in those 
values is contained within the width of the curve. 














